#
# load_and_merge_xprs.R
# Load and combine an arbitrary number of eXpress "results.xprs" outputs
# and generate a cumulative bundle_id taking into account all of the bundled
# targets in each file.
# 20140610
#
# Copyright (C) 2014 Shawn Driscoll, sdriscoll@salk.edu
#
# This program is free software; you can redistribute it and/or
# modify it under the terms of the GNU General Public License
# as published by the Free Software Foundation; either version 2
# of the License, or (at your option) any later version.
# 
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
# 
# You should have received a copy of the GNU General Public License
# along with this program; if not, write to the Free Software
# Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
#
# Description:
# Merges several eXpress 'xprs' files and creates a final overall 
# bundle id for the de-novo bundling of targets. By default a connection
# between any two targets only needs to occur once (i.e. in a single xprs file)
# however it might be useful to set the min.occ value to the minimum number of 
# replicates in any one condition or to just simply choose an arbitraty number
# of samples that makes sense at the time.
#
# DEPENDS ON igraph package. install with:
# install.packages("igraph")
#
require(igraph)
load_and_merge_xprs <- function(flist, min.occ=1, max.bundle.size=100) {
	
	# helper function
	explode <- function(s, delim=";") {
		unlist(strsplit(s, delim))
	}

	# load the data from the provided file list
	lload <- lapply(flist, function(f) {
				message(sprintf("loading %s", f))
				d <- read.delim(f, stringsAsFactors=F)
				rownames(d) <- d[, 2]
				return(d)
			})
	
	# build graph of the target names linked by bundle
	
	# vertex names will be all of the target ids
	message("generating target id vector")
	vi <- unique(unlist(lapply(lload, function(a) a[, 2])))
	
	message("finding all bundled targets")

	# split all target ids by bundle id in each file
	tmp <- lapply(lload, function(a) split(a[, 2], a[, 1]))
	# concatenate
	tmp2 <- do.call(c, tmp)
	rm(tmp)
	# reduce to only those with more than one target name
	idx <- which(sapply(tmp2, length) > 1)
	tmp2 <- tmp2[idx]

	# remove those with length greater than maximum bundle size. apparantely eXpress
	# sometimes bundles massive numbers of transcripts together. maybe due to 
	# excessive mismatch allowance during alignment.
	sizes <- sapply(tmp2, length)
	idx <- which(sizes > max.bundle.size)
	if(length(idx) > 0) {
		message(sprintf("dropping %d bundles larger than %d. max bundle was %d", 
			length(idx), max.bundle.size, max(sizes)))
		tmp2 <- tmp2[-idx]
	}

	message("building graph")

	# build matrix of edges
	tmp <- lapply(tmp2, function(a) {
		combn(a, 2)
	})
	tmp <- t(do.call(cbind, tmp))

	# reduce the size of this by summing together redundant edges (this won't catch those in
	# opposite orders but those will reduce later in the graph)
	id <- paste(tmp[, 1], tmp[, 2], sep=":")
	id.t <- tapply(id, id, length)
	da.id <- data.frame(t(sapply(names(id.t), explode, delim=":")), id.t)
	names(da.id) <- c("from", "to", "weight")
	# make graph
	g <- graph.data.frame(da.id, vertices=vi, directed=F)
	
	# simplify the graph down to combine redundant edges
	g <- simplify(g, edge.attr.comb=sum)
	# drop edges below minimum occurence threshold
	cix <- which(E(g)$weight < min.occ)
	if(length(cix) > 0) {
		g <- g-E(g)[cix]
	}

	message("building output")

	# create new bundle ids from the graph's connected component clusters
	gi <- clusters(g)$membership
	
	# make a bundle table
	bundles <- data.frame(bundle_id=gi, target_id=V(g)$name, stringsAsFactors=F)
	bundles <- bundles[order(by=bundles[, 1]), ]
	
	# build matrices of different columns from each sample
	cnames <- names(lload[[1]])[3:ncol(lload[[1]])]
	ltables <- vector(mode="list", length=length(cnames))
	names(ltables) <- cnames
	for(n in cnames) {
		ltables[[n]] <- sapply(lload, function(a) a[bundles$target_id, n])
		rownames(ltables[[n]]) <- bundles$target_id
	}
	
	lout <- list(bundles=bundles, flist=flist)
	lout <- c(lout, ltables)
	
	return(lout)
}
